# !/usr/bin/python3
# -*- coding: utf-8 -*-
# @Time : 2021/9/1 10:16
# @Author : Yu LiXinQian
# @Email : yulixinqian805@gmail.com
# @File : PCA_TE.py
# @Project : My_MLF


import numpy as np
from time import perf_counter
from scipy.io import loadmat
from scipy.stats import chi2, f
import matplotlib.pyplot as plt

print("="*80)
t1 = perf_counter()
print("Waiting loading Data!")
TE_data = loadmat("Data/TE_Data.mat")
d00_te = TE_data['d00_te']  # 960*52
d21_te = TE_data['d01_te']
X00_1 = d00_te[:, 0:22]
X00_2 = d00_te[:, 41:52]
X21_1 = d21_te[:, 0:22]
X21_2 = d21_te[:, 41:52]
X00 = np.hstack((X00_1, X00_2))
X21 = np.hstack((X21_1, X21_2))
t2 = perf_counter()
print("Successfully Executed in {:.2f}ms!".format((t2-t1)*1000))

NumSamples, NumVariables = d00_te.shape

print("There are the data for d00_te, d21_te!")
print("There are {} samples, and {} variables per sample!".format(NumSamples, NumVariables))
print("="*80)

t3 = perf_counter()
print("Normalize the Data!")
x_mean = X00.mean(axis=0, keepdims=True)
x_std = X00.std(axis=0, keepdims=True)

X00 = X00 - X00.mean(axis=0, keepdims=True)
X00 = X00 / X00.std(axis=0, keepdims=True)

X21 = X21 - x_mean
X21 = X21 / x_std

t4 = perf_counter()
print("Done in {:.2f}ms!".format((t4-t3)*1000))
print("="*80)

t5 = perf_counter()
print("Singular value Decomposition!")
U, S, VT = np.linalg.svd(X00/np.sqrt(NumSamples-1))
V = VT.T
SigmasMat = np.diag(S)
t6 = perf_counter()

print("U:{}".format(U.shape))
print("SigmasMat:{}".format(SigmasMat.shape))
print("V:{}".format(V.shape))
print("Done in {:.2f}ms!".format((t6-t5)*1000))
print("="*80)

t7 = perf_counter()
print("Computing the contribution of eigenvalues!")
eigenvalues = np.power(S, 2)
per_contribution = (eigenvalues / np.sum(eigenvalues)).reshape(33, 1)
sum_contribution = 0
for i in range(len(per_contribution)):
    sum_contribution += per_contribution[i, 0]
    if sum_contribution >= 0.8:
        PCs = i+1
        break
t8 = perf_counter()
# print("per_contribution:{}".format((np.around(per_contribution.flatten()*100, 1))[0:PCs]))
print("No.{} 's contribution is over 80%!, and it's {:.2f}%!".format(PCs, sum_contribution*100))
print("per_contribution:{}".format((np.around(per_contribution.flatten()*100, 1))[0:PCs]))
print("Done in {:.2f}ms!".format((t8-t7)*1000))
print("="*80)

t9 = perf_counter()
P = V[:, 0:13]
Score = np.dot(X00, P)
X00_recover = np.dot(Score, P.T)
E = X00 - X00_recover
SPE_Normal = np.sum(np.power(E, 2), axis=1, keepdims=True)
g = SPE_Normal.var(axis=0)/(2*SPE_Normal.mean(axis=0))
h = 2*np.power(SPE_Normal.mean(axis=0), 2)/SPE_Normal.var(axis=0)
SPE_limit = g*chi2.ppf(0.99, h)
t10 = perf_counter()
print("The control value of SPE_limit:{:.2f}".format(SPE_limit[0]))
print("Done in {:.2f}ms!".format((t10-t9)*1000))
print("="*80)

t11 = perf_counter()
T2_Normal = np.zeros((NumSamples, 1))
for i in range(NumSamples):
    for j in range(PCs):
        T2_Normal[i, 0] += np.power(Score[i, j], 2)/ eigenvalues[j]
T2_limit = PCs*(NumSamples-1)*(NumSamples+1)/NumSamples/(NumSamples-PCs)*f.ppf(0.99, PCs, NumSamples)
t12 = perf_counter()
print("The control value of T2_limit:{:.2f}".format(T2_limit))
print("Done in {:.2f}ms!".format((t12-t11)*1000))
print("="*80)

t13 = perf_counter()
Score21 = np.dot(X21, P)
X21_recover = np.dot(Score21, P.T)
E21 = X21 - X21_recover
SPE_21 = np.sum(np.power(E21, 2), axis=1, keepdims=True)
T2_21 = np.zeros((NumSamples, 1))
for i in range(NumSamples):
    for j in range(PCs):
        T2_21[i, 0] += np.power(Score21[i, j], 2)/ eigenvalues[j]

t14 = perf_counter()
print("="*80)

print("Make Plots!")
t15 = perf_counter()

plt.figure(figsize=(16,10))
plt.subplot(2, 2, 1)
plt.plot(T2_Normal)
plt.plot(np.repeat(T2_limit, NumSamples), color='red')
plt.title("T2 Of Normal Data")
plt.subplot(2, 2, 2)
plt.plot(SPE_Normal)
plt.plot(np.repeat(SPE_limit, NumSamples), color='red')
plt.title("SPE Of Normal Data")

plt.subplot(2, 2, 3)
plt.plot(T2_21)
plt.plot(np.repeat(T2_limit, NumSamples), color='red')
plt.title("T2 Of Normal Data")
plt.subplot(2, 2, 4)
plt.plot(SPE_21)
plt.plot(np.repeat(SPE_limit, NumSamples), color='red')
plt.title("SPE Of Normal Data")
t16 = perf_counter()
print("Done in {:.2f}ms!".format((t16-t15)*1000))
plt.savefig("Output/PCA_TE/test.png")
plt.show()


